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Quantum Monte Carlo Calculations of Dihydrogen Binding Energetics on Ca Cations: an 
Assessment of Errors in Density Functionals for Weakly Bonded Systems 
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We investigate the binding of single and quadruple hydrogen molecules on a positively charged Ca ion. 
By comparing with benchmark quantum Monte Carlo (QMC) calculations we demonstrate wide variability in 
other more approximate electronic structure methods including common density functionals. Single determi- 
nant QMC calculations find no binding at short range by approximately 0.1 eV for the quadruple hydrogen 
molecule case, for a fixed hydrogen bond length of 0.77 Angstrom. Density functional calculations using com- 
mon functionals such a LDA and B3LYP differ substantially from the QMC binding curve. We show that use 
of full Hartree-Fock exchange and PBE correlation (HFX+PBEC) obtains close agreement with the QMC re- 
sults, both qualitatively and quantitatively. These results both motivate the use and development of improved 
functionals and indicate that caution is required applying electronic structure methods to weakly bound systems 
such as hydrogen storage materials based on metal ion decorated nanostructures. 
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PACS numbers: 

Our ability to accurately predict molecular adsorption ener- 
gies is of widespread importance in the physical, chemical, 
and materials sciences. Technologically, the adsorption of 
small molecules on semiconductors or metals is an essential 
step in many catalytic or energy storage related areas. In the 
case of hydrogen storage, the strength of the adsorption can 
determine the suitability a material for practical application: 
if the binding is too high, release of the hydrogen will be dif- 
ficult at moderate operating temperatures, while if the binding 
is too weak, storage of the hydrogen will be ineffective. 

Calculating the energetics of hydrogen adsorption is a dif- 
ficult task and requires highly accurate quantum mechanics 
based calculations. If the structure and eventually the dynam- 
ics of the adsorption process are to be accurately modeled, 
the potential energy surface of the adsorbant and adsorbate 
must be accurately simulated over a length scale of at least 
5 Angstrom. Density functional theory (DFT) based meth- 
ods are the most widely applied electronic structure methods 
for studies of hydrogen storage materials. However, in prac- 
tice the DFTs are not only approximate, but are also rarely 
benchmarked in the non-bonding and weakly bonding config- 
urations vital for hydrogen storage. 

Motivated in part by recent discussions and discrepancies 
for DFT predictions of hydrogen adsorption on alkaline earth 
metals jjJ43|], we have performed extensive Quantum Monte 
Carlo (QMQIJ] calculations for H2 (dihydrogen) adsorption 
on the Ca 1+ system. The ion's charge models the scenario 
where the ion is absorbed on graphene[5, 6]. QMC provides 
an accurate and unbiased reference to compare against ap- 
proximate but more computationally affordable approaches. 
We concentrate on the interaction of hydrogens with a single 
ion, as opposed to a system with a substrate, since the geome- 
tries are easily and unambiguously specified and the systems 
are already sufficient to demonstrate substantial differences in 
predicted binding energies and overall shape of the binding 



The interaction of one or many H2 molecules with Ca 1+ 
will undoubtedly involve several effects: charge transfer, po- 
larization, and potential long range dispersion (or van der 
Waals) interactions. To accurately model these systems, first 
principles calculations should be able to accurately account 
for all these effects with little reliance on, e.g. error cancela- 
tion. For example, van der Waals interactions interactions are 
naturally and accurately included within QMC appro aches 17(- 
llOll . but are absent from common DFTs. 

In the following we (i) describe our QMC methodology, (ii) 
present benchmark results for the cases of single H2 adsorp- 
tion on Ca 1+ , and (iii) since in actual use scenarios additional 
hydrogen molecules will be present we also results for quadru- 
ple H2 adsorption. These systems are constructed identically 
to those of Ref. (JLfl . Finally, (iv) we summarize our findings. 

Quantum Monte Carlo — The QMC method allows for a 
very efficient and accurate solution of the Schrodinger equa- 
tion. In contrast with many electronic structure methods, 
QMC methods involve only well controlled approximations. 
Although their computational prefactor is often large, for 
small and medium sized molecular systems energetics close to 
chemical accuracy can be obtained, e.g. Refs. IUIII12II . These 
properties make QMC methods ideal for benchmark studies 
and for the assessment of computationally cheaper but more 
approximate methods. 

QMC methods are wavefunction based and the most im- 
portant input is the trial many body wavefunction. In Varia- 
tional Monte Carlo (VMC) a direct variational evaluation of 
the energy of a trial wavefunction is performed using impor- 
tance sampled Monte Carlo integration. VMC calculations 
therefore suffer from a potentially very strong dependence 
on the input wave function and any prior assumptions about 
the electronic structure, but have the advantage that the actual 
many-body wavefunction is obtained and can be analyzed. In 
fixed-node diffusion quantum Monte Carlo (DMC), the low- 
est energy state consistent with the zeros (nodes) of the trial 



wavefunction is projected. This projection greatly reduces the 
dependence of the final energy on the input trial wavefunction 
compared to VMC. In practice, very accurate results are ob- 
tained byDMC for a wide variety of molecular and solid state 
svstemsi4l7l 4ld,ll3ill4ll . Due to the increased robustness we 
concentrate on DMC results in this study. 

For our purposes, the only significant approximations in 
DMC calculations are (i) the use of pseudopotentials and (ii) 
the fixed-node approximation and consequent dependence on 
the nodal surface of the input trial wavefunction. The first 
approximation introduces systematic errors via the approx- 
imate treatment of core-valence interactions and via the lo- 
cality approximation [L5J necessary to evaluate the non-local 
pseudopotentials in DMC. We minimize these errors by using 
a small Ne core for the Ca pseudopotential 1 1611 and very high 
quality trial wavefunctions. We use the same pseudopotentials 
in all our calculations to ensure a fair comparison between all 
methods: the same Hamiltonian is solved in our QMC, DFT, 
and quantum chemical calculations. 

To minimize the nodal errors in our DMC calculations we 
also use multideterminant trial wavefunctions obtained from 
configuration interaction calculations that are subsequently re- 
optimized via the energy minimization method ! 1 ill . This ap- 
proach is a significant advance over conventional applications 
of DMC where the nodal surface of the trial wavefunction 
consists of only a single Slater determinant determined by a 
less accurate theory such as DFT: the nodal errors are sys- 
tematically reduced to near chemical accuracy ll ill 12ll when 
sufficient statistics can be obtained and the multideterminant 
expansion is large enough. Previous studies have shown that 
(i) for light molecules single determinant DMC yields results 
similar in accuracy CCSD(T) with the aug-cc-pVQZ basis 
set lU3ll . (ii) these errors are further reduced with multideter- 
minant methods, e.g. Illll - ll3ll , and (iii) pseudopotential er- 
rors are small and less significant than the nodal error in these 
systemsClH- 

In principle, modern trial wavefunction optimization 
methods 1 1 ill 1211 can produce DMC results nearly independent 
of the input provided sufficiently flexible trial wavefunction 
forms are adopted. Here we validate our single determinant 
nodal surface results using large configuration interaction ex- 
pansions of many determinants. 

In the following calculations we use trial wavefunctions 
consisting of a weighted sum of Slater determinants multi- 
plied by a two-body Jastrow factor. The Slater determinants 
consist of orbitals determined by GAMES S i 1711 DFT or com- 
plete active space self consistent field (CASSCF) calculations 
expanded in the large ANO-VTZ gaussian basis set llol . The 
two-body Jastrow factor does not change the nodal surface but 
acts to enforce the electron-electron cusp condition, greatly 
improving the overall quality of the trial wavefunctions. For 
the QMC calculations we used the QWALK code|QJl. Mul- 
tideterminant QMC calculations used up to 370 determinants, 
where we took all CASSCF determinant of squared magni- 
tude greater than 0.01. Energy minimization was performed 
starting from the truncated CASSCF results. We note that the 



DMC energies always lie substantially below the pure quan- 
tum chemistry results. Pseudopotentials were derived in the 
soft Hartree-Fock formalism 1 1611 . An average DMC popula- 
tion of ~ 30000 walkers and a small time step of 0.005 a.u. 
was used. The largest DMC calculations used 0(1000) pro- 
cessor hours per energy point. 

Results for single H2 absorption — Figure [T] shows our cal- 
culated binding energy curve for the hydrogen dimer on Ca 1+ , 
with the molecular bond oriented perpendicularly to the line 
of approach. Since the calculation of forces is not well devel- 
oped in QMC, for each distance from the Ca 1+ ion we com- 
puted energies for all methods with a fixed bond length of 0.77 
Angstrom, corresponding to the value found near binding in 
quantum chemical calculations |1]. Only small changes in the 
fully relaxed value are seen over all distances, indicating that 
the trends in the binding at fixed bond length are representa- 
tive of the relaxed case. For computational simplicity, we treat 
the energy of the system at z = 4.6 Angstrom as representing 
fully separated unbound system. 

Our results show that while the potential energy surface 
varies quantitatively between the methods, for a single dimer 
the general trends given are qualitatively similar for many of 
the methods, with a single minimum. However, unrestricted 
Hartree-Fock (UHF) and second order M0ller-Plesset per- 
turbation theory (MP2) calculations show negligible binding. 
The DMC data shows a minimum around z = 3.1 Angstrom, 
and binding of ~ 0.025 eV. 

Comparing the density functional results against the DMC 
energy curve we find that B3LYP ll9ll functional gives a rela- 
tively good agreement, with minimum at z =< 2.9 Angstrom. 
However, the LSDA functional [20] significantly overbinds by 
at least 0.1 eV, while the PBE functional 1 2 ill lies midway be- 
tween the B3LYP and LSDA values. None of the functionals 
results in false energetic minima in the binding energy curve, 
however the distance of the minimum energy varies by 0.5 
Angstrom over these functionals. Calculations using Hartree- 
Fock exchange combined with PBE correlation (HFX+PBEC) 
(similar to Refs. |22|, |23| except with 100% exchange) very 



closely resemble the DMC results; analysis and possible rea- 
sons for the apparent accuracy are discussed after the four di- 
hydrogen results. 

Results for four H2 absorption — Figure|2]shows our calcu- 
lated binding energy curve for four hydrogen dimers on Ca 1+ . 
In this system the hydrogens are pinned in a planar geometry, 
ninety degrees apart in D4 symmetry, with molecular bonds 
oriented perpendicularly to the line of approach (inset in Fig- 
ure O. 

The binding energies obtained with single determinant 
DMC display a minimum around 2.2 Angstrom. However, 
comparing these energies with those over 3 Angstrom clearly 
shows the minimum to be a local metastable minimum: there 
is no binding of four H2 molecules in this planar geometry 
at short range at the single determinant DMC level. To test 
the accuracy of these calculations we also used multideter- 
minant wavefunctions determinants, initially obtained by re- 
stricted active space RAS(9,37) calculations. The binding is 
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FIG. 1: Calculated binding energy for a single hydrogen molecule 
(H2) approaching Ca 1+ . The molecule to oriented with the bond 
perpendicular to the line of approach (see inset). Results are shown 
for the unrestricted Hartree-Fock (UHF), the local spin density ap- 
proximation (LSDA), PBE and B3LYP density functionals, second 
order M0ller-Plesset perturbation theory (MP2). We also show Dif- 
fusion Quantum Monte Carlo (DMC) and density functional results 
calculated using exact exchange combined with PBE correlation 
(HFX+PBEC). The DMC calculations use a single determinant of 
B3LYP orbitals and a Jastrow factor for the trial wavefunction and 
nodal surface. The hydrogen molecule bond length is held fixed at 
0.77 Angstrom. The lines are a guide to the eye. Error bars are 
smaller than the size of the DMC symbols. 



shifted to higher energies by ~ 0.1 eV indicating that the sin- 
gle determinant results and nodal surface are robust. 

DFT based calculations show clear energy minima around 
the 2.2 Angstrom distance indicating significant binding of 
the H 2 molecules. Although the depth of binding varies, sim- 
ilar behavior is obtained for LSDA, PBE and B3LYP. The 
same energetic ordering is observed as for the single hydro- 
gen case, with LSDA displaying greatest binding. By con- 
trast, UHF calculations display only a slight minimum around 
2.3 Angstrom. We also include quantum chemical results 
from RAS and complete active space (CAS) calculations in 
Fig|2] The larger active space calculations reduce the calcu- 
lated binding energy, moving the quantum chemical results 
towards the DMC results. Perturbative theory results (MP2) 
also show no overall binding. A clear transition between states 
of Ai g and B2 g symmetry is observedyj between 2.5 and 3.0 
Angstrom depending on the underlying theory. 

Our qualitative conclusion of no binding for the four hy- 
drogen molecule case with fixed 0.77 Angstrom bond length 
is in qualitative agreement with previous quantum chemical 
results[ 1]. However our more extensive basis sets and more 
rigorous DMC calculations reveal that that the energy scale for 
any potential binding is very small, only a few tenths of an eV. 
Strikingly, even allowing a generous estimate of 0.2eV resid- 
ual systematic errors in our DMC calculations, any eventual 
binding will remain small (order O.leV) whereas the LSDA 
and PBE functionals predict binding energies one order of 



magnitude larger. 

We also tested energies obtained from density functional 
theory using using Hartree-Fock exchange combined with 
PBE correlation, gradually increasing the fraction of ex- 
change. As for the single hydrogen molecule case, the cal- 
culated binding curve accurately follows the DMC data over 
all distances. However, only a 100% contribution fully repro- 
duced the DMC data; lesser contributions smoothly interpo- 
lating between the DMC and PBE results. Given the wide 
variation seen for other functionals, this is noteworthy, and 
also indicates the primary source of error in the other den- 
sity functional predictions: We argue, in accordance with Ref. 
12411 . that the semi-local functionals provide good description 
of static correlation and exchange in this system. On the other 
hand, dynamic or long-range exchange is mostly absent in 
these functionals. Since our system seems not to have a strong 
multi-reference character (based on our RAS and CAS calcu- 
lations) the dynamical part of exchange must play a domi- 
nant role. This last part is well described only in HF or in 
exact-exchange functionals such as optimized effective poten- 
tial (OEP) method. Therefore including full dynamical part of 
a exchange together with the correlation from LSDA, PBE, or 
B3LYP gives a good description. This observation is mostly 
independent from the type of the correlation used. 

Summary — We have performed benchmark quantum 
Monte Carlo calculations of hydrogen molecule binding on 
the Ca cation. Density functional calculations vary widely in 
the predicted binding energy. Density functional calculations 
using Hartree-Fock exchange well reproduce the Monte Carlo 
results, suggesting a route - if not universal - to predictive and 
accurate calculations in this and related systems. We hope that 
our results will help further motivate the development of im- 
proved functionals. As metal ions have been proposed as hy- 
drogen binding centers in new hydrogen storage materials, we 
strongly recommend caution in applying density functional 
methods to these systems. Appropriate benchmarking using 
quantum chemical or quantum Monte Carlo techniques is re- 
quired. 
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factor for the trial wavefunction and nodal surface (black triangles). At 2.2 Angstrom separation we also compute the binding energy in DMC 
using Energy Minimized Restricted Active Space (RAS) multideterminant wavefunctions. At this distance we also show quantum chemical 
results for RAS(9,37) (red triangle), RAS (9,31) (blue circle), and complete active space CAS(9,18) (orange circle). HFX+PBEC indicates 
density functional results calculated using exact exchange combined with PBE correlation. The hydrogen molecule bond lengths are held fixed 
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